# ------------------ importing simulation library ---------------------- #
from actinNetwork import *
from numpy import arccos,sqrt,median,argmax,std,unique,savetxt,vstack
from numpy.random import normal,choice
from scipy.special import erf
from scipy.stats import zscore
from matplotlib.pyplot import subplots,bar,xticks,axhline,setp
import csv

# display plots in notebook
%matplotlib inline

# ------------------------------- parameters for simulation ----------------------------- #

# initial number of filaments N, and monomer size dN, and range in nm
N = 200; dN = 5.5; filRange = 1000.0;

# running time T, precision dt, and frame rate ds
T = 10.0; dt = 0.001; ds = 2.0; 

# velocity and width of polymerisation profile 
Fo = 0.2; dw = 5*dN
vo = 250.0

# ------------------------------- velocity pertubation ---------------------------------- #

# time and duration of perturbation
tau = 4.0; dtau = 0.1

# velocity change during perturbation
dF = 0.2

# velocity perturbation as heaviside function centered at tau
def F(t) :
    return Fo - dF * heavisidePi( t-tau-sqrt(dtau), 2.0*sqrt(dtau) )


# ------------------------------- rate assumptions ---------------------------------- #

# define polymerisation rate in terms of
# critical angle given unperturbed velocity
phiCrit = 50.0 * pi/180

# rates for travelling heavisides
lambdaRate = vo/(dN*cos(phiCrit))
betaRate = 0.5
kappaRate = 1000.0

# defining rate functions
def rLambda(x,y,theta,t) :
    return lambdaRate*heavisidePi(y,dw)
   
def rBeta(x,y,theta,t) :
    return betaRate*heavisidePi(y,dw)

def rKappa(x,y,theta,t) :
    return kappaRate*heavisideTheta(-y-dw/2)+1.0*heavisidePi(y,dw)

plot([t for t in arange(0,T,dt)],[F(t) for t in arange(0,T,dt)])
plot([t for t in arange(0,T,dt)],[Fo + exp(-(t-1.065*tau)**2) for t in arange(0,T,dt)])

[<matplotlib.lines.Line2D at 0x7fd1bcb3aa50>]

# binning cross x direction
xEdges1 = linspace(0,2000,27)
xCenters1 = (xEdges1[:-1] + xEdges1[1:]) / 2.0

xEdges2 = linspace(0,2000,14)
xCenters2 = (xEdges2[:-1] + xEdges2[1:]) / 2.0

# binning across phi direction
phiEdges1 = array([0,20,50,90])*pi/180
phiCenters1 = ( phiEdges1[1:]+phiEdges1[:-1] ) / 2.0

phiEdges2 = linspace(-pi/2,pi/2,10)
phiCenters2 = ( phiEdges2[1:]+phiEdges2[:-1] ) / 2.0

# function that calculates relevant statistics given a network object that has already been run
def getStatistics(networkObject,xEdges1,xEdges2,phiEdges1,phiEdges2) :
    
    # getting positions
    xFil = networkObject.getPositions( networkObject.Monomers )
    xBranch = networkObject.getPositions( networkObject.Branches )
    xCap = networkObject.getPositions( networkObject.Caps )

    # angular orientation of filaments
    phiFil = networkObject.getAngles ( networkObject.Monomers )
    phiBranch = networkObject.getAngles ( networkObject.Branches )
    phiCap = networkObject.getAngles ( networkObject.Caps )

    # calculating normalised filament denisty in x direction
    filamentDensity = histogram( xFil.T[1], bins = xEdges1)[0]
    filamentDensity = filamentDensity/float(max(filamentDensity))

    # calculating zscore branching and capping rates in x direction
    branchRate = zscore( histogram( xBranch.T[1], bins = xEdges1)[0] )
    capRate = zscore( histogram( xCap.T[1], bins = xEdges1)[0] )

    # binning monomers into x-slices
    filMasks = [ ( x <= xFil.T[1] ) * ( xFil.T[1] <= y ) for x,y in zip(xEdges2[:-1],xEdges2[1:]) ]

    # and using these bins to calculate angule histograms at each slice
    filNumber = array([ histogram( abs(unique(phiFil[mask])), bins = phiEdges1)[0] for mask in filMasks ]).astype(float)
    phiDensity = array([ histogram( phiFil[mask], bins=phiEdges2)[0] for mask in filMasks ])

    # calculation of order parameter
    orderEdges = array([-50,-30,-10,10,30,50])*pi/180
    orderHist = [ histogram(phiFil[mask],bins=orderEdges)[0] for mask in filMasks ]
    orderParam = [ ( N[2]-(N[0]+N[-1])/2.0 ) / float( N[2]+(N[0]+N[-1])/2.0 ) for N in orderHist ]
    
    # return relevant statistics
    return filamentDensity,branchRate,capRate,filNumber,phiDensity,orderParam

#xSteady = vstack([n.getPositions(n.Frontier).T[0],uniform(low=0.0, high=dw/2.0, size=len(n.Frontier))]).T;
#dxSteady = array([ [dN*sin(theta),dN*cos(theta)] for theta in n.getAngles(n.Frontier) ])

# ------------------------------- initialisation ---------------------------------- #

# initial distribution of network
try : # using previous positions and angles if possible
    
    # getting previous positions and angles
    #xInit = n.getPositions(n.Frontier) - array(len(n.Frontier)*[[0.0,p(n.tElapsed)]])
    #dxInit = array([ [dN*sin(theta),dN*cos(theta)] for theta in n.getAngles(n.Frontier) ])
    
    xInit = xSteady
    dxInit = dxSteady
    
except : # otherwise initialise uniformly
    
    # sampling from uniform distribution given N
    thetaDist = uniform(-pi,pi,size=N)
    xDist = uniform(low=0.0, high=filRange, size=N)
    yDist = uniform(low=0.0, high=dw/2.0, size=N)
    
    # passing samples to inital variables
    xInit = array([ [x,y] for x,y in zip(xDist,yDist) ])
    dxInit = array([ [dN*sin(theta),dN*cos(theta)] for theta in thetaDist ])
    
# initialise network object
n = network(rLambda, rBeta, rKappa,
            xSeed = xInit, dxSeed = dxInit,
            branchSigma = 15.0*pi/180.0,
            forceDirection = True, recordHistory = True)

# ------------------------------- simulation ---------------------------------- #

# evolve network up to specified time
#n.exportData(dt,ds,n.tElapsed+T,Fext=1.0e-1)
#n.evolve(dt,n.tElapsed+T,Fext=1.0)

# while networks grows and up until time tFinal
while n.nBarbed != 0 and n.tElapsed <= T :

    # evolve
    n.timeStep(dt, Fext = F(n.tElapsed) )

### ------------------------------- Filament Plot ------------------------------ ###

# getting positions
xFil = n.getPositions( n.Monomers )
xBranch = n.getPositions( n.Branches )
xCap = n.getPositions( n.Caps )

figure(figsize=(20,5))

# plotting all monomers as different colour points 
plot(xFil.T[1],xFil.T[0],'g',marker=".",linewidth=0,ms=2)
plot(xBranch.T[1],xBranch.T[0],'#2737ff',marker=".",linewidth=0,ms=4)
plot(xCap.T[1],xCap.T[0],'#ff0000',marker=".",linewidth=0,ms=4)

# plotting options
xlabel(r"Distance, $x$ / nm",fontsize=18)
ylabel(r"Distance, $y$ / nm",fontsize=18)
ylim(0,n.xBoundary); xlim(0,2000);

# ------------------------------- calculate statistics ---------------------------------- #

filamentDensity,branchRate,capRate,filNumber,phiDist,orderParam = getStatistics(n,xEdges1,xEdges2,phiEdges1,phiEdges2)

### ------------------------------- Density Plot ------------------------------ ###

f,ax1 = subplots()
ax2 = ax1.twinx()

f.set_figheight(5)
f.set_figwidth(20)

# plotting rates and density
ax1.plot( xCenters1, filamentDensity, '#00ff00',mec='#00ff00',marker="o",mew="4",ms=12,linewidth=4)
ax2.plot( xCenters1, branchRate, '#2737ff',mec='#2737ff',marker="o",mew="4",ms=12,linewidth=4)
ax2.plot( xCenters1, capRate, '#ff0000',mec='#ff0000',marker="o",mew="4",ms=12,linewidth=4)

# plotting options
ax1.set_xlabel(r"Distance, $x$ / nm",fontsize=16)
ax1.set_ylabel(r"Percentage Filament Denisty, $\rho(x)$",fontsize=16)
ax2.set_ylabel(r"Standardised Rate, $\omega(x)$", fontsize=16)
ax1.set_ylim(0,1); ax1.set_xlim(0,)

### ------------------------------- Capping Ratios ------------------------------ ###

figure(figsize=(20,5))

# plot capping ratios
plot( xCenters2, filNumber.T[0],'k',
     mec='k',marker="o",mew="4",ms=12,linewidth=4,
     label=r"0$^{\circ}$- 20$^{\circ}$")
plot( xCenters2, 2.0/3.0*filNumber.T[1],'orange',
     mec='orange',marker="o",mew="4",ms=12,linewidth=4,
     label=r"20$^{\circ}$- 50$^{\circ}$")
plot( xCenters2, 1.0/2.0*filNumber.T[2],'0.7',
     mec='0.7', marker="o",mew="4",ms=12,linewidth=4,
     label=r"50$^{\circ}$- 90$^{\circ}$")

# plotting options
xlabel(r"Distance, $x$ / nm",fontsize=16); l = legend(title=r"Angle $\phi$",fontsize=16)
ylabel(r"Filament Number $N(x\vert\phi)$",fontsize=16)
setp(l.get_title(),fontsize=16)
xlim(0,)

### ------------------------------- Angle Order ------------------------------ ###

f,ax1 = subplots()
ax2 = ax1.twinx()

f.set_figheight(5)
f.set_figwidth(20)

# plotting angle distributions
for x,y in zip([ linspace(x-40,x+40,len(phiCenters2)) for x in xCenters2 ],phiDist) :
    ax1.bar( x, y, width=10, color="gray", edgecolor="k", linewidth=2, align="center")

# and order parameter
ax2.plot( xCenters2, orderParam, "r.-", marker="o", mec="r",mew="4",ms=12,linewidth=4, label="Order Parameter")
ax2.axhline( 0, color="k", linewidth=1);

# plotting options
ax1.set_ylabel(r"Filament Density, $\rho(\phi\vert x)$",fontsize=16);
ax2.set_ylabel(r"Order Parameter,$\alpha(x)$",fontsize=16);
ax1.set_xlabel(r"Distance, $x$ / nm",fontsize=16)
xlim(0,);

# initialising arrays to gather results
FD = []
BR = []
CR = []
AD = []
OP = []

# number of times to run the simulation
nRuns = 20

# looping over nRuns
for run in range(nRuns) :
    
    # ------------------------------- initialisation ---------------------------------- #

    # initial distribution of network
    try : # using previous positions and angles if possible

        # getting previous positions and angles
        # xInit = n.getPositions(n.Frontier) - array(len(n.Frontier)*[[0.0,p(n.tElapsed)]])
        # dxInit = array([ [dN*sin(theta),dN*cos(theta)] for theta in n.getAngles(n.Frontier) ])
        
        #xSteady = vstack([n.getPositions(n.Frontier).T[0],uniform(low=0.0, high=dw/2.0, size=len(n.Frontier))]).T;
        #dxSteady = array([ [dN*sin(theta),dN*cos(theta)] for theta in n.getAngles(n.Frontier) ])
        
        xInit = xSteady
        dxInit = dxSteady

    except : # otherwise initialise uniformly

        # sampling from uniform distribution given N
        thetaDist = uniform(-pi,pi,size=N)
        xDist = uniform(low=0.0, high=filRange, size=N)
        yDist = uniform(low=0.0, high=dw/2.0, size=N)

        # passing samples to inital variables
        xInit = array([ [x,y] for x,y in zip(xDist,yDist) ])
        dxInit = array([ [dN*sin(theta),dN*cos(theta)] for theta in thetaDist ])

    # initialise network object
    n = network(rLambda, rBeta, rKappa,
                xSeed = xInit, dxSeed = dxInit,
                branchSigma = 15.0*pi/180.0,
                forceDirection = True, recordHistory = True)

    # ------------------------------- simulation ---------------------------------- #

    # while networks grows and up until time tFinal
    while n.nBarbed != 0 and n.tElapsed <= T :

        # evolve
        n.timeStep(dt, Fext = F(n.tElapsed) )

    ### ----------------------------- Processing Data --------------------------- ###

    filamentDensity,branchRate,capRate,filNumber,phiDist,orderParam = getStatistics(n,xEdges1,xEdges2,phiEdges1,phiEdges2)

    FD += [filamentDensity]
    BR += [branchRate]
    CR += [capRate]
    AD += [filNumber]
    OP += [orderParam]

### ------------------------------- Filament Plot ------------------------------ ###

# getting positions
xFil = n.getPositions( n.Monomers )
xBranch = n.getPositions( n.Branches )
xCap = n.getPositions( n.Caps )

figure(figsize=(20,5))

# plotting all monomers as different colour points 
plot(xFil.T[1],xFil.T[0],'g',marker=".",linewidth=0,ms=5,alpha=0.2)
plot(xBranch.T[1],xBranch.T[0],'#2737ff',marker=".",linewidth=0,ms=6,alpha=0.5)
plot(xCap.T[1],xCap.T[0],'#ff0000',marker=".",linewidth=0,ms=6,alpha=0.5)

# plotting options
xlabel(r"Distance, $x$ / nm",fontsize=18)
ylabel(r"Distance, $y$ / nm",fontsize=18)
ylim(0,n.xBoundary); xlim(0,2000);

### ------------------------------- Density Plot ------------------------------ ###

f,ax1 = subplots()
ax2 = ax1.twinx()

f.set_figheight(5)
f.set_figwidth(20)

# plotting rates and density

ax1.errorbar( xCenters1, mean(array(FD),axis=0),yerr=std(array(FD),axis=0), color='#00ff00',mec='#00ff00',marker="o",mew="4",ms=0,linewidth=4)
ax2.errorbar( xCenters1, mean(array(BR),axis=0),yerr=std(array(BR),axis=0), color='#2737ff',mec='#2737ff',marker="o",mew="4",ms=0,linewidth=4)
ax2.errorbar( xCenters1, mean(array(CR),axis=0),yerr=std(array(CR),axis=0), color='#ff0000',mec='#ff0000',marker="o",mew="4",ms=0,linewidth=4)

# plotting options
ax1.set_xlabel(r"Distance, $x$ / nm",fontsize=18)
ax1.set_ylabel(r"Percentage Filament Denisty, $\rho(x)$",fontsize=18)
ax2.set_ylabel(r"Standardised Rate, $\omega(x)$", fontsize=16)
ax1.set_ylim(0,1); ax2.set_ylim(-5,5); ax1.set_xlim(0,2000);

f,ax = subplots()
f.set_figheight(5)
f.set_figwidth(20)

# plot capping ratios
ax.errorbar( xCenters2, mean(array(AD),axis=0).T[0],yerr=std(array(AD),axis=0).T[0],color='k',
     mec='k',marker="o",mew="4",ms=0,linewidth=4,
     label=r"0$^{\circ}$- 20$^{\circ}$")
ax.errorbar( xCenters2, 2.0/3.0*mean(array(AD),axis=0).T[1],yerr=2.0/3.0*std(array(AD),axis=0).T[1],color='orange',
     mec='orange',marker="o",mew="4",ms=0,linewidth=4,
     label=r"20$^{\circ}$- 50$^{\circ}$")
ax.errorbar( xCenters2, 1.0/2.0*mean(array(AD),axis=0).T[2],yerr=1.0/2.0*std(array(AD),axis=0).T[2],color='0.7',
     mec='0.7', marker="o",mew="4",ms=0,linewidth=4,
     label=r"50$^{\circ}$- 90$^{\circ}$")

# plotting options
xlabel(r"Distance, $x$ / nm",fontsize=18); legend(fontsize=16)
ylabel(r"Percentage Angular Denisty $\rho(x\vert\phi)$",fontsize=18)
xlim(0,2000);

f,ax1 = subplots()
ax2 = ax1.twinx()

f.set_figheight(5)
f.set_figwidth(20)

# plotting angle distributions
for x,y in zip([ linspace(x-80,x+80,len(phiCenters2)) for x in xCenters2 ],phiDist) :
    ax1.bar( x, y, width=20, color="none", edgecolor="k", linewidth=2, align="center")

# and order parameter
ax2.errorbar( xCenters2, mean(array(OP),axis=0),yerr=std(array(OP),axis=0), color="r", linewidth=4, label="Order Parameter")
ax2.axhline( 0, color="k", linewidth=1);

# plotting options
ax1.set_ylabel(r"Number of Monomers, $N(\phi\vert x)$",fontsize=16);
ax2.set_ylabel(r"Order Parameter, $\alpha(x)$",fontsize=16);
ax1.set_xlabel(r"Distance, $x$ / nm",fontsize=16)
xlim(0,2000);

# ---------------------------------------save everything!! ----------------------- #
savetxt("X-FD-BR-CR-correlated.csv",
array([xCenters1,
mean(array(FD),axis=0),std(array(FD),axis=0),
mean(array(BR),axis=0),std(array(BR),axis=0),
mean(array(CR),axis=0),std(array(CR),axis=0) ]).T,delimiter=",")

savetxt("X-AD-correlated.csv",array([xCenters2,
1.0/1.0*mean(array(AD),axis=0).T[0],1.0/1.0*std(array(AD),axis=0).T[0],
2.0/3.0*mean(array(AD),axis=0).T[1],2.0/3.0*std(array(AD),axis=0).T[1],
1.0/2.0*mean(array(AD),axis=0).T[2],1.0/2.0*std(array(AD),axis=0).T[2] ]).T,delimiter=",")

savetxt("X-OP-correlated.csv",array([xCenters2,
mean(array(OP),axis=0),std(array(OP),axis=0) ]).T,delimiter=",")

# binning monomers into x-slices
phiFil = n.getAngles ( n.Monomers )
xFil = n.getPositions( n.Monomers )
filMasks = [ ( x <= xFil.T[1] ) * ( xFil.T[1] <= y ) for x,y in zip(xEdges2[:-1],xEdges2[1:]) ]

b = open("SAD-correlated.csv", 'w')
a = csv.writer(b)
data = [[str(x*180.0/pi) for x in phiFil[mask]] for mask in filMasks ]
a.writerows(data)
b.close()

for i,x in enumerate(OP[:2]) :
    
    plot(x,color=str(i/float(len(OP))))
    
axhline( 0, color="k", linewidth=1)

<matplotlib.lines.Line2D at 0x7fa770caa450>

### ------------------------------- Filament Plot ------------------------------ ###

# getting positions
xFil = n.getPositions( n.Monomers )
xBranch = n.getPositions( n.Branches )
xCap = n.getPositions( n.Caps )

figure(figsize=(80,20))

# plotting all monomers as different colour points 
plot(xFil.T[1],xFil.T[0],'g',marker=".",linewidth=0,ms=12,alpha=0.2)
plot(xBranch.T[1],xBranch.T[0],'#2737ff',marker=".",linewidth=0,ms=15,alpha=0.5)
plot(xCap.T[1],xCap.T[0],'#ff0000',marker=".",linewidth=0,ms=15,alpha=0.5)

# plotting options
axis("off");

n.nBarbed

438